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A thermal ablation model for silicates is proposed. The model includes the mass losses through the 
balance between evaporation and condensation, and through the moving molten layer driven by surface 
shear force and pressure gradient. This model can be applied in ablation simulations of the meteoroid or 
glassy Thermal Protection Systems for spacecraft. Time-dependent axi-symmetric computations are 
performed by coupling the fluid dynamics code, Data-Parallel Line Relaxation program, with the material 
response code, Two-dimensional Implicit Thermal Ablation simulation program, to predict the mass lost 
rates and shape change. For model validation, the surface recession of fused amorphous quartz rod is 
computed, and the recession predictions reasonably agree with available data. The present parametric 
studies for two groups of meteoroid earth entry conditions indicate that the mass loss through moving 


molten layer is negligibly small for heat-flux conditions at around 1 MW/cm’. 


Nomenclature 


= normalized ablation rate 

= specific heat, J/kg-K 

= heat transfer coefficient, J/m?-s-K 
= mass transfer coefficient, kg/m?-s 
= diffusion coefficient, m/s 

= total internal energy, J/m° 

= enthalpy, J/kg 

= recovery enthalpy, J/kg 


= thermal conductivity of solid, W/m-K 
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= molecular weight, kg/kmol 

= mass flux, kg/m?-s 

= mass loss rate of moving molten layer, kg/m?-s 
= mass loss rate due to evaporation, kg/m7?-s 

= pressure, Pa 

= saturated vapor pressure, Pa 

= heat flux, W/m? 

= in-depth conductive heat flux at surface, W/m? 
= convective heat flux at surface, W/m? 

= radiative heat flux at surface, W/m? 

= translation-vibration energy exchange rate, W/m? 
= radius, cm 

= corner radius of the model, cm 

= nose radius of the model, cm 


= universal gas constant, 8.3143 J/ mol-K 


= temperature, K 


= maximum allowed change on surface temperature, K 
= environment temperature, K 


= flow velocity, m/s 

= speed of moving grid, m/s 

= thermal velocity, m/s 

= Cartesian coordinate system, m 

= source term of gas-surface interactions, kmol/m?-s 
= absorptance 

= evaporation/condensation sticking coefficient 


= effective thickness of moving liquid layer, cm 


American Institute of Aeronautics and Astronautics 


E = emissivity 


p = total density, kg/m° 

X = blowing reduction coefficient 
u = viscosity, Pa-s 

Oo = Stefan-Boltzmann constant, W/(m?-K*) 
T = shear force, Pa 

subscripts 

] = unblown 

C = char 

l, J = direction components 

= molten layer 

n = surface species 

S = gas species 

V = vapor or vibrational energy 

Ww = wall 


I. Introduction 


Chondrites are the most common type of meteoroid that falls to Earth. Most chondrites are rich in silicate materials. 
To determine the physical and chemical nature of the meteoroid, one must understand the interaction between the 
meteoroid and the atmosphere. In the typical meteor ablation model [1], ablation begins when the surface of the 
meteoroid reaches the boiling temperature, and then the temperature is assumed to remain constant. The classic meteor 
ablation theory assumes that mass loss is proportional to the kinetic energy imparted to the meteoroid. Meteoroids 
may lose mass not only though evaporation, but also through spraying of the molten layer on the surface. To simulate 
this effect, the classic ablation model includes a mass loss term which is an inverse tangent function of the difference 
between surface temperature and melting point. Another commonly used mass loss model for meteoroid is the mass 
loss rate computed using Langmuir evaporation assuming the condensation rate 1s negligibly small. [2] 

In this work, a generalized thermal ablation model for silicate materials is proposed. This model includes the mass 


losses through the balance between evaporation and condensation, and through the moving of molten layer driven by 
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surface shear force and pressure gradient. This model is currently implemented in the Two-dimensional Implicit 
Thermal AblatioN simulation program, TITAN [3], which was developed for the computation of Thermal Protection 
System materials under hypersonic entry conditions. The ablation characteristics of silicates are the same as those of 
silica and metal, which are commonly used materials for the Thermal Protection Systems for NASA’s space missions. 
Hence, this generalized ablation model can be applied for the meteoroid, or any glassy and metallic components of 
spacecraft Thermal Protection Systems. The time dependent axi-symmetric simulations for mass loss and shape 
change predictions are performed by coupling the fluid dynamics code, Data-Parallel Line Relaxation (DPLR) [4], 
and the material response code, TITAN. The computational predictions are compared with available data for model 
validation.[13] Parametric studies under meteoroid entry conditions are performed to gain further understanding on 


the ablation characteristics of chondritic materials in order to provide guidance for future work. 


II. Thermal Ablation Model 


The mechanism of thermal ablation for silicate materials includes mass loss rate due to the motion of molten layer 
driven by surface shear force and surface pressure gradient, m,, and mass loss rate due to the balance between 
evaporation and condensation at the surface, m,. 

When subjected to aerothermodynamics heating and shear forces, the molten material will form a flowing liquid 
layer over the solid surface. The equations governing the liquid layer are nearly the same as for the case of an 
incompressible gas boundary layer, the difference being a variation of viscosity with temperature. For an axi- 
symmetric geometry, if the inertia force in momentum equations can be neglected, the momentum and continuity 


equations of molten layer can be integrated to obtain the mass loss rate due to the motion of liquid layer [6]: 


rin, = 24 {~ (7,6? - 2¥59)} (1) 


r ax “Uw ax 


The viscosity of molten layer usually increases exponentially with a decrease in temperature 1.e. 


a 


p=er? (2) 
where a and b are constants. The “effective thickness” of a moving liquid layer, 0, is defined as the distance at which 
the viscosity increases by one exponential factor and it’s assumed that the thickness of the moving liquid layer is thin 
in comparison with the thickness of thermal layer. 

The evaporation rate of a surface species can be computed from its saturated vapor pressure, Pf, and its 


condensation rate from the gas phase can be computed if the species partial surface pressure, P, , is known [7] : 
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thy = Yin YM, (PE — Ph) (3) 
where V,, is the thermal speed of surface species n, and y is the evaporation/condensation (or sticking) coefficient. 

If the ambient gas pressure is very low, the condensation rate is thus negligible. Then the computation of m,, can 
be simplified to 

ty = in YMyPe (4) 

The saturated vapor pressure can be computed using a chemical equilibrium code, such as MAGMA [8]. The 
MAGMA code was developed specifically for the composition of meteoroids or earth crust. Forsterite, Mg2SiO4(s), 
is the most common silicate material which can be found in chondrite meteorites. In MAGMA, Mg2SiQO.(1) is assumed 
to undergo pseudo-reaction to form pseudo-species, MgQO(1) and SiO2(1) which then becomes vapor, MgO(g) and 
S102(g), as the temperature increases. The computed saturated vapor pressures as function of surface temperature for 
MgO (1) and SiQ2(1) are presented in Fig. 1. 

If a continuous flow is formed over the surface, the species partial pressure at the surface needs to be calculated 
using a finite-rate reacting flow solver. The gas phase chemical species and their chemical reactions for ablation 
products have to be considered in the reacting flow simulation [9]. A finite-rate reacting flow simulation with all the 
chemical species of ablations products is complicated and extremely CPU-time consuming. However, if the chemical 
equilibrium condition can be established at the surface, the net mass loss rate due to evaporation/condensation process 
can be pre-calculated, i.e., the B’ tables, using a chemical equilibrium code such as MAT [10]. The computed Bz of 
silica, forsterite, and a hypothetical chondrite in air as a function of surface temperature and pressure are presented in 
Figs. 2. The hypothetical chondrite is assumed to have the elemental mass fractions of O = 0.4465, Si = 0.2018, Mg 
= 0.1613, S = 0.0225, and Fe = 0.1679, based on the chemical composition of an ordinary chondrite. The computed 
Bc curves show that the boiling temperatures of chondrite at various surface pressures are more similar to those of 
forsterite than silica. Figure 3 presents the mole fractions of the most abundant twelve gas species at the surface of 
chrodrite as a function of B¢ at a surface pressure of 100 atm, assuming that the solid surface is in chemical equilibrium 
with the adjacent air. For the meteor earth entry conditions that we are interested in this work, the maximum surface 
heat flux reaches around the order of 1 MW/cm/?, and the maximum surface pressure is 100 atm. A significant amount 
of evaporation (i.e. a high Bz value) is thus expected, and the most abundant ablation product should be SiO followed 
by O2 and MgO. 


The mass loss due to evaporation/condensation for a chemical equilibrium surface is written as 
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m, = CyBe (5) 
The total mass loss rate due to thermal ablation, m,, is thus expressed as 


Me = My “f= Mm, (6) 
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Figure 1: Saturated vapor pressures for SiO2 and MgO. 
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Figure 3: Mole fractions of the most abundant twelve gas 
species at the chondrite surface with the surface pressure of 
100 atm, if the surface is in chemical equilibrium with the 
adjacent air. 





III. Computational Methods 


In this section the relevant equations of the TITAN and DPLR codes are described and the coupling procedure 
employed to obtain the solutions is presented. 
A. Fluid-Solid Coupling 
Coupling between the material response and flow codes is required for many multidimensional ablation 
simulations, because the magnitudes and distributions of the surface heat flux, pressure, and shear force are very 
sensitive to shape change. The integration between TITAN and DPLR is based on a loosely coupled approach 
following the flowchart presented in Fig. 4 [11]. In this approach, all CFD solutions are computed for an unblown 
surface but use the ablated shape. Chemical equilibrium surface thermochemical interactions and blowing reduction 
effects are incorporated in the material-response code by use of ablation tables, the surface energy balance with heat 


transfer coefficient, and the blowing reduction parameter. 
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The initial flow field and its associated surface heat flux, pressure, and shear are first computed by DPLR. For 
each surface point, heat transfer coefficient, Ca, shear force, T,,,, and pressure, p, are calculated. Then shear force, 
pressure, and heat transfer coefficient are passed as boundary conditions to TITAN. A time-accurate ablation and 
thermal response computation is performed by TITAN. When the maximum surface temperature or recession meets a 
pre-specified limit, TITAN stops its computation and outputs the location of the ablated surface. A new CFD grid is 
then generated based on this ablated surface and estimated shock location, and a new steady-state flow solution is then 
calculated by DPLR. The surface pressure, shear force, and heat transfer coefficients are calculated for this new 
solution. The updated values for pressure, shear force, and heat transfer coefficients are input to TITAN for another 
run of time-accurate ablation and in-depth thermal response. This procedure is repeated until TITAN reaches the 
specified final time. In each CFD run, the outer boundary of grid is aligned with the shock using an internal subroutine 
in DPLR, and the cell Reynolds number equal to | is enforced. The TITAN-DPLR coupling procedure is written in a 


UNIX/Linux script file. The test cases presented in this paper are computed on a Linux cluster system. 


TITAN 


(time marching from 


t* to t*+ At,) , | ) DPLR 


(steady-state) 





Figure 4. DPLR/TITAN shape change coupling methodology. 


B. Material Ablation and Thermal Response Code 
The governing equations in TITAN, which include energy conservation, a three-component decomposition model, 


and the surface energy balance, are solved with a moving grid system to simulate the shape change due to surface 
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recession [3]. The equations are converted from a Cartesian coordinate system (x,y) into a general body-fitted 
coordinate system (é, 7) and then discretized using a finite-volume method. A moving grid technique was 
implemented to solve simulations with surface recession. 

The internal energy balance for a silicate material is a transient thermal conduction equation without pyrolysis 
terms: 

Plp <= V- (KVT) + pepv VT (7) 

The individual terms in Equation (7) may be interpreted as follows: the rate of storage of sensible energy, the net rate 
of thermal conductive heat flux, and the convection rate of sensible energy due to coordinate system movement. 

Conditions at the ablating surface are determined by the aerothermal environment and by chemical interactions 
between the boundary layer gas, the ablation products, and the chemical constituents of the surface material. TITAN 
employs a convective transfer-coefficient form of the surface energy balance for silicate: 

Cy (Hy, — Ay) + 1c(he — hw) + AwOrw — FEw(Tw — Te) — Icona = 0 (8) 

Here PeueCy is the convective heat transfer coefficient, H; is the recovery enthalpy, and all other quantities are defined 
at the ablating surface. The first term in Equation (8) is the convective heat flux; the second term represents the 
chemical energy released (or absorbed) by ablation; the third and fourth terms are radiation absorption and emission, 
respectively; and the final term is the rate of heat conduction into the TPS. This equation has been simplified from the 
general form by assuming equal diffusion coefficients within the boundary layer and equal Stanton numbers for heat 
and mass transfer. These are standard assumptions for most entry environments and TPS materials. If these 
assumptions are not applicable, a more general form of the surface energy balance may be used. It is also assumed 
that flow-field radiation is reflected or absorbed at the surface, but not transmitted. 

A blowing correction accounts for the reduction in heat transfer coefficient due to the injection of gases from 
pyrolysis and surface ablation into the boundary layer. The blowing correction equation used by TITAN for silicate 
iS: 


Cy _ In(1+24B')  ,, _ me 


Cie er ee (9) 


Here / is the blowing reduction parameter, and Cy/Cm 1s the ratio of the blown (ablating) to the unblown (non- 
ablating) heat transfer coefficients. For laminar flow, 4 is 0.5 or higher, depending on the geometry and the ratio of 


molecular weights of the injected and boundary-layer-edge gas. For transitional or turbulent flow, smaller values of / 
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are used. B’ is the non-dimensional mass blowing rate at the surface. Unless noted otherwise, the laminar blowing 
reduction parameter of 0.5 is used for the calculations presented herein. By assuming that the momentum boundary 
layer thickness is the same as the thermal boundary layer thickness, the blockage of surface shear force for an ablating 
surface is proportional to that of heat transfer coefficient, and is approximated by the following equation from 


Reference [6]: 


i aa, (OT (10) 


Tw1 CH1 


C. Fluid Dynamics Code 

The Navier-Stokes solver, DPLR, is used to estimate the hypersonic aero-heating distribution over a blunt body 
[4]. The governing equations may be characterized as representing a flow-field in thermal and chemical non- 
equilibrium. The DPLR code solves the time-dependent conservation equations of mass, momentum, and energy 


within the flow-field. The species mass conservation equation is given by: 


Ops 0 


0 
De. + ax; pst) = ~ 9q, (sPsi) + Ws (11) 


The total momentum conservation is written as: 





0 0 es OTij 
=, (Pui) + a (puju;) = a (12) 
The vibrational and total energy equations are written as: 
dEy , a ee 
at t ax, Boy) = — 5, ey) + Cr (13) 
dE. a a a a 
oe bx ((E + p)u;) =-— ax (q; + vj) ae (u;ti;) = eey ax, Dsilts (14) 


In DPLR, the Gauss-Seidel line relaxation method is modified to enable fast convergence for viscous flow with 


high parallel efficiency on massively parallel computers. 


IV. Results 


The computations presented in this section contain two parts. In the first part, validation of the proposed ablation 
model for silicate materials is performed. Silica is considered a silicate, although it is the special case with no negative 


charge and no need for counter-ions. The ablation characteristics of silica are assumed to be the same as those of other 
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silicate materials. Since the thermal properties of quartz are well studied, quartz is one of the best candidate materials 
for our ablation model validation. The ablation data for quartz published by Brogan [13] is used to evaluate the 
prediction accuracy of proposed ablation model for silicates. In the second part, the parametric studies of ablation 
characteristics for a hypothetical chondritic material under two groups of meteoroid entry conditions are presented. 
The focus of parametric studies is to explore how the various meteoroid entry environments [12] affect the 
contributions of moving molten layer and evaporation to the total surface ablation. 

Part I: Model Validation — Quartz ablation 

The test model used for validation is a 0.88-cm radius fused amorphous quartz rod. Two types of quartz materials 
are investigated, a clear quartz and an opaque quartz.[5][13] The stagnation point pressure is 0.13 atm, and the initial 
cold wall heat flux for a hemisphere with the same nose radius is 950 W/cm’. The thermal properties of silica suggested 
in Reference 5 are adapted in the computation. The ablating surface is assumed to be in chemical equilibrium with the 
adjacent arc stream. Thus the surface evaporation rate is computed by using the pre-calculated ablation B; tables 
shown in Fig. 2a. The rate of moving molten layer is computed by the DPLR/TITAN coupled simulation. The total 
recession rate is the sum of evaporation and moving molten layer rates as described in Eqs. 1, 5, and 6. 

The comparison of stagnation point recession vs time between computation and measurement is presented in 
Fig.5. The solid lines stand for the predictions and the symbols represent data. The total recession of clear quartz 
model is about twice as much as that of opaque quartz model at the end of 25 s. Also note that for the first 15 s, the 
predicted recession for both clear and opaque quartz models are higher than data. This is probably because the 
assumption of chemical equilibrium surface is inconsistent with the low surface catalytic efficiency of glassy materials 
when the surface mass blowing rate is relatively low. After around 15 s, the agreement between computation and 
measurement is generally better as compared to that in the first 15 s. The predicted recession rates due to evaporation 
(vapor) and moving molten layer (liquid) for clear and opaque models are presented in Figs. 6a and 6b, respectively. 
The evaporation rate is a strong function of surface temperature. Thus, the vapor recession rate reaches its steady state 
condition in about 5 s. The rate of moving liquid depends on the depth of thermal layer. Therefore, before the thermal 
layer can grow and establish a steady-state condition, the vapor recession rate is dominant. As soon as the thermal 
layer establishes a steady thickness at around 12-13 s, the liquid moving rate becomes higher than the vapor recession 
rate for both clear and opaque quartz models. Fig. 7a shows the predicted surface temperature history for both models. 


It takes about 5 s for the surface temperature to reach a steady state condition. Hence the predicted surface evaporation 
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rates shown in Figs 6 also take about the same amount of time to reach a steady state condition. The predicted steady- 
state surface temperature is around the boiling temperature (~2700 K) of silica at pressure of 0.13 atm. The surface 
temperature of opaque model is about 50 K below that of the clear quartz model because of the difference on surface 
emissivity. Figures 7b and 7c show in-depth temperature histories along the center-line of the model at 5, 10, 15, 20, 
and 25 s. The profiles of in-depth temperatures at 15, 20, and 25 s are similar for both clear and opaque models. This 
implies that the depth of thermal penetration has established a steady condition before 15 s, resulting from the balance 
between surface recession rate and thermal diffusion rate. This trend is consistent with the predicted moving liquid 


rates shown in Figs. 6. 
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Figure 5: Comparisons of computed and measured | Figure 6a: Predicted stagnation point recession rates for 
stagnation point recession for clear and opaque quartz | the clear quartz model. 
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Part II: Parametric Studies 

In this section, we study two groups of conditions for meteoroid atmospheric entry [12]. The main purpose is to 
investigate how the variation of entry conditions affect the contributions due to moving molten layer and evaporation 
to the total surface recession. The primary goal is to guide the direction of future work and to identify the relevant 
material properties which have to be obtained to improve the prediction accuracy. The first group is with an initial 
nose radius of 5 m, a short total entry time of 5 s, a heat flux at order of 1 MW/cm/”, and a surface pressure at 100 atm. 
The second group is with a smaller initial nose radius of 5 cm, a total entry time of 30 sec, a moderate maximum heat 
flux at order of 1 kW/cm/, and a surface pressure at 1 atm. A triangular heat pulse is assumed in these simulations. 
Since a detailed thermal response model for a chondritic material is not currently available, the thermal properties 
measured at room temperature of 300 K are used in these parametric studies. The density is assumed to be 3000 kg/m’, 
the thermal conductivity is 4.5 W/m-K, the specific heat is 690 J/kg-K, and the surface emissivity is 0.6. Clearly, these 
assumptions may have some negative impact on the accuracy of predicted surface recession and in-depth thermal 
response. However, they should be adequate for performing these parametric studies. The ablating surface is assumed 
to be at chemical equilibrium condition because of high surface pressure, and thus the ablation B; tables shown in 
Fig. 2c and Eq. 5 are used in the computation of surface evaporation rate. The rate of moving liquid is calculated 
through the coupling between DPLR and TITAN. This hypothetical chondrite is assumed to be a piece of 
homogeneous solid material. Spallation and fragmentation are not considered in our ablation simulation. 

In the first group of computation, the surface heating boundary condition is a 5-sec triangular pulse. The maximum 
total heat fluxes of 1, 1.5, and 2 MW/cm/ , and the same maximum surface pressure of 100 atm are imposed in the 
boundary conditions. This is to examine the sensitivity of heat flux on the predicted surface recession rate. Figure 8a 
presents the predicted stagnation point evaporation rates and moving liquid rates for these three triangular heat pulses. 
The dashed lines are the evaporation rates and the solid lines are the moving liquid rates. The predicted evaporation 
rates are directly proportional to the imposed surface heat fluxes. The peak surface evaporation rate for the maximum 
heat flux of 2 MW/cm? reaches about 35 cm/s, and that for the 1 MW/cm/? is about 50% of that for the 2 MW/cm’. 
The predicted moving liquid rates appear to be negligibly small as compared with the evaporation rates. The predicted 
total surface recession (dashed lines) and the surface temperature (solid lines) are shown in Fig. 8b. The difference on 


predicted surface temperature among three heating conditions 1s relatively small. This is because the imposed surface 
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heat fluxes for these three conditions are high enough for the surface to reach the boiling temperature of around 4500 
K at the surface pressure of 100 atm. Hence, the total amount of material removed through evaporation 1s directly 
proportional to the total amount of energy available at the surface. The in-depth temperature profiles along the central 
line for 2 MW/cm* condition at times equal to 1, 2,3, 4, and 5 s are shown in Fig. 8c. At this heating rate, the vaporized 
surface is receding at a speed far much higher than the thermal diffusion speed. Thus, an extremely sharp temperature 
drop at the front surface is predicted through the entire 5-second period of exposure. Therefore, the interior temperature 
is unaffected by the surface heating and remains almost unchanged during the entire entry. Figure 8d is the in-depth 
temperature profile at time equal to 4 s. It is plotted in a very fine length scale to resolve the detail structure of this 
very sharp temperature drop. It indicates that the depth of thermal penetration is only about 0.1 mm at time equal to 4 
s. This is why the rate of moving liquid layer is insignificant under these extreme heating environments. 

Figures 9 present the recession predictions for three different maximum surface pressures (100, 200, and 300 atm) 
at the same maximum heat flux of 1 MW/cm7. The predicted evaporation rates are not sensitive to the variation of 
surface pressure. The pressure change just slightly alter the evaporation rate, and in comparison the rate of moving 
liquid layer remains negligible. However, the variation of surface pressure does affect the surface temperatures, which 
are shown in Fig. 9b. The boiling temperature increases as surface pressure increases, and thus higher surface pressure 
results in higher surface temperature. The condition with a higher surface temperature/pressure requires extra small 
amount of energy to reach its boiling temperature, and hence the high pressure condition has slightly lower recession 
than the low pressure condition does. At entry conditions similar to what studied in this group, surface radiation and 
in-depth conduction (the last two terms in Eq.8) are much smaller than the total surface heating. 

The total energy required for evaporation is the sum of the energy required to reach the boiling point from the 
initial temperature and the enthalpy difference between liquid and vapor at the boiling point. Since the rate of moving 
liquid can be ignored, the total mass loss can be estimated from the total surface heat load if the total energy required 
for evaporation is known. Hence, in order to calculate the total recession for these extremely high heating conditions 
the following material properties are required: the specific heat of solid/liquid as a function of temperature and the 
heat of formation at the reference temperature. The rest of the thermal properties are considered to be not as important 


in these recession calculations. 
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In the second group of study, a 30-sec triangular heat pulse is imposed as the heating boundary condition. The 
maximum heat fluxes of 1, 1.5, and 2 kW/cm? with the same maximum surface pressure of | atm are used in the 
computation. Figure 10a presents the comparison of predicted stagnation point evaporation rates (dashed lines) and 
moving molten layer rates (solid lines) among these three heat fluxes. As expected, these predicted recession rates are 
far smaller than those presented in the first group, and are dominated by the moving liquid. The shape of predicted 
surface evaporation rate is again similar to the shape of the imposed triangular heat pulse. This is because the 
evaporation rate is mostly sensitive to the surface temperature. The predicted maximum evaporation rate for the 2 
kW/cm/ heat pulse is just about 0.1 mm/s. Since the surface heat fluxes are not high enough to obtain a large scale 
evaporation as in the first group and the depth of thermal penetration is thick enough to form a significant molten 
layer, the moving molten layer becomes the major contributor to the surface recession. The evaporation process 
contributes less than 5% of total surface recession. The predicted total surface recession (dashed lines) and the surface 
temperature (solid lines) are shown in Fig. 10b. The peak surface temperatures reach some where in-between 3000 to 
3200 K. The boiling temperature is around 3200 K at the surface pressure of | atm. The predicted total surface 
recessions are not directly proportional to the total heat load, because majority of the surface recession is the result of 
moving liquid layer. The in-depth temperature profiles along the axis of symmetry for 2 kW/cm? condition at times 
equal to 5, 10, 15, 20, 25, and 30 s are shown in Fig.10c. The predictions indicate that the depth of thermal penetration 
is around 0.5 cm at time equal to 5 s, and grows to about | cm at time equal to 10 s. Thus, a liquid layer forms at the 
front surface where the temperature is above the melting point, and it moves when the surface shear force or pressure 
gradient is significant. 

Figures 11 present the predicted recession rates at the maximum heat flux of 2 kW/cm/’, including the sensitivity 
of recession to the shear force. The rate of moving liquid layer is dependent on the magnitude of the surface shear 
force. Therefore, doubling the shear increases the moving liquid rate, and halving the shear decreases the moving 
liquid rate. On the other hand, the evaporation rate is a function of surface temperature, and it is not directly affected 
by the surface shear value. However, the surface temperature changes slightly when the moving rate of molten layer 
changes as indicated in Fig. 11b. Consequently, the evaporation rate is also affected by the surface shear, but this 


perturbation is far smaller than the liquid moving rate. 
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Variation of surface pressure can affect the moving liquid rate and evaporation rate. This is because the boiling 
point depends on the surface pressure, and the moving liquid layer is driven by the derivative of the surface pressure 
as described in Eq. 1. Figures 12a and 12b present the predictions of surface recession rates under three different 
maximum surface pressures of 0.5, 1, and 1.5 atm, and at the same maximum heat flux of 2 kW/cm/?. The predictions 
shown in Figs 12 are consistent with the theory. Decreasing or increasing the surface pressure changes the boiling 
point thus affecting both the predicted evaporation rate and the predicted surface temperature. Since the moving rate 
of liquid layer is proportional to the surface pressure gradient, scaling down or up the surface pressure decreases or 
increases the liquid moving rate. 

Unlike the conditions in the first group, the in-depth penetration of thermal diffusion is significant and the 
formation of molten layer is the dominant contributor to surface recession in the second group of conditions. Therefore, 
to accurately predict the recession rate for conditions similar to those in the second group, all the thermal properties, 
such as thermal conductivity, specific heat, and emissivity as functions of temperature for chondritic material, are 
equally important and have to be obtained. Additionally, the viscosity of liquid layer as a function of temperature has 
to be measured for the modeling of the motion of liquid layer. A simple heat of evaporation or ablation approach 


discussed in the first group would not be appropriate for the second group. 
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V. Conclusions 

The formulation of a thermal ablation model for silicate materials is presented and discussed. The recession 
computation is conducted through the coupling between the Navier-Stokes solver, DPLR, and the material simulation 
code, TITAN. Model validation has been accomplished by using the recession data for a 0.88-cm radius fused 
amorphous quartz rod. The agreement on total surface recession between computation and measurement is generally 
good. Parametric studies for two groups of meteoroid entry conditions have been performed. The results show that 
moving molten layer and in-depth thermal penetration are negligibly small for the extremely high heat flux (order of 
1 MW/cm/’) entry conditions, and a simple Heat of Ablation Model should be adequate for the recession prediction. 
Additionally, for conditions with heat flux around the order of 1 kW/cm/?, the surface recession is mostly driven by 
moving molten layer, and thus a detailed thermal response model of chondrite and the viscosity of its molten layer are 
required for the ablation simulation. In order to further refine the present ablation model for silicates, it is proposed 
that new arcjet tests should be performed to obtain more data points. These arcjet tests would be under conditions 


where the moving molten layer is the dominant ablation mechanism. 
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